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ABSTRACT 


Acoustic  tomography  signal  transmissions  in  the  Monterey  Bay  is 
modelled  using  the  time-domain  parabolic  equation  method  of  Collins  and 
Westwood  (1991).  Comparison  of  the  model  output  with  the  measured 
arrival  structures  obtained  in  Monterey  Bay  in  1988  shows  that  this 
Fourier  synthesis  can  produce  good  agreement  with  data.  Furthermore, 
identification  of  the  measured  modal  arrivals  is  possible  by  decomposing 
the  PE  model  output  into  individual  normal  modes.  Modal  identification 
provides  for  the  application  of  tomography  in  shallow  water. 
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L  INTRODUCTION 


In  December,  1988,  phase-encoded  acoustic  signals  were  transmitted 
continuously  for  four  days  in  the  waters  off  Monterey,  California  to  investigate 
the  feasibility  of  using  tomography  for  ocean  interior  and  surface  monitoring  in  a 
coastal  environment.*'^  Sound  arrival  structures  were  measured  from  a  source 
(center  frequency  of  224  Hz  and  bandwidth  of  16  Hz),  located  in  the  deep,  881  m, 
sound  channel  to  several  receivers  on  the  shallow,  100  m,  adjacent  continental 
shelf  (Fig,  1).  The  path  from  the  source  to  the  receiver  at  Station  J  is  shown  in  Fig. 
2.  Every  16  received  1.9736  second  time  periods  were  coherently  averaged.  The 
squared  magnitude  of  every  other  averaged  sequence  is  displayed  as  a  function 
of  time  in  Fig.  3,  for  Station  J  starting  on  December  14.  Stable  arrival  structures 
were  observed  over  the  duration  of  the  experiment  for  this  path, 

A  successful  application  of  tomography  depends  upon  the  stability, 
resolvability  and  identifiability  of  the  acoustic  arrivals.  If  acoustic  ray  or  mode 
arrival  structures  are  stable,  and  the  individual  arrivals  have  sufficient  signal  to 
noise  ratio  (SNR)  and  are  resolved  in  time,  then  the  travel  times  of  the  rays  or 
modes  can  be  estimated  over  the  duration  of  the  experiment.  It  is  the  nature  of  the 
acoustic  tomography  technique  that  travel  time  perturbation  time  series  be 
calculated  as  the  differences  between  the  measured  arrival  times  and  modeled  or 
predicted  times  of  arrivals.  It  is  these  time  series  which  are  the  data  for  the  inverse 
problem.  Therefore,  to  calculate  the  perturbation,  the  identification  problem  in 
tomography  must  be  solved  where  actual  acoustic  ray  or  mode  arrivals  are 
associated  with  modeled  or  predicted  arrivals.  The  arrivals  in  deep  water  can  be 
identified  using  ray  theory.^'^ 
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Ray  tracing  techniques  had  mixed  success  in  modeling  the  propagation  in 
the  shallow  water  of  Monterey  Canyon^,  and  in  the  Florida  Straits.*  In  both  these 
cases,  the  number  of  rays  present  in  shallow  water  became  large  enough  to 
prevent  the  identification  of  individual  rays  using  a  single  receiver.  In  addition. 
Smith  found  that  acoustic  rays  were  extremely  sensitive  to  the  exact  bathymetry 
specification  of  the  Monterey  Canyon  continental  slope.^ 

A  full  wave  approach  provides  a  more  realistic  model  of  the  transmitted 
signal  in  the  presence  of  large  bathymetric  changes.  The  parabolic  equation 
method  developed  by  Collins^’^^  was  chosen  for  this  purpose  due  its  energy 
conservation  scheme,  and  its  ability  to  efficiently  incorporate  higher  order  terms 
in  its  approximation  of  the  Helmholtz  equation.  The  steep  continental  slope  gives 
rise  to  higher  order  terms  of  significant  magnitude,  and  Collins  showed  that 
energy  conservation  is  essential  when  modeling  slopes  in  excess  of  ten  degrees.^ 
Sound  arrival  identification  is  accomplished  by  the  decomposition  of  the 
pressure  field  generated  by  the  PE  model  into  normal  modes.  Additionally,  this 
process  provides  increased  physical  understanding  of  mode  coupling  and  modal 
sound  transmission.  Tomographic  inverse  methods  using  normal  mode  travel  times 
have  been  formulated  within  the  context  discussed  by  Munk  and  Wunsch  for  the 
ocean  interior,  and  Miller  et  al.  for  the  surface. 
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Figure  1.  Location  of  Source  (A)  and  Receivers  (B  through  L-2)  for  the 
Monterey  Bay  Tomography  Experiment. 


Figure  2.  Three«diinensional  section  of  Monterey  Bay  with  source  (A) 

Station  J  locations  marked. 
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Figure  3.  Station  J  measured  1.9736  second  segment  averages  vs  time 
received  (decimal  hours  starting  on  December  14, 1988  PST).  Every  other  16 
segment  average  is  plotted  for  visual  clarity. 
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n.  PRESSURE  CALCULATION  USING  FOURIER  SYNTHESIS 

1.  A  T  ?  ME  DOMAIN  APPROACH 

The  receiver  time  domain  pressure  signal  can  be  determined,  under  the 
assumption  of  linearity  and  time  invariance,  as  a  convolution  of  the  source 
pressure  and  the  impulse  response  of  the  ocean  waveguide;^'* 


rr,z,t)  =  | 


p(ro ,  zo ,  t-x)h(ro ,  zq  ,  Tr ,  z,  x)dT 


where  p  is  pressure,  ro  is  the  source  range,  rr  is  the  receiver  range,  zq  is  source 
depth,  Zr  is  receiver  depth  and  t  is  time. 

The  Fourier  transform  of  the  convolution  integral  in  the  time  domain  is  a 
multiplication  in  the  frequency  domain: 

P(rr,z,f)  =  P(ro,zo,f)*  H(ro,zo,rr,z,f) 


Where  P  and  H  are  the  Fourier  transforms  of  p  and  h  respectively  such  that: 


p(rr,z,  t)  = 


roo 

P(rr, 

J-OO 


z, 


6 


h(ro,Zo,rr,z, 


t)  =  f  H(ro,  zo,  Tr ,  z,  f)e-i2"fkif 


Therefore,  the  pressure  in  the  time  domain  can  also  be  determined  from  the  inverse 
Fourier  transform  of  (2): 

foo 

p(Xr ,  Z,  t)  =  P(ro  ,Zo ,  f)H(ro ,  Zq  ,  ,  z,  f) 

J-oo 


2.  DETERMINATION  OF  THE  TRANSFER  FUNCTION 

The  Helmholtz  equation  governs  the  sound  pressure  disturbances  excited 
by  a  point  harmonic  source.  Therefore,  the  transfer  function  can  be  obtained  by 
solving  the  Helmholtz  equation  at  different  frequencies.  Allowing  for  density 

variation,  the  Helmholtz  equation  is: 

pV-  (-^Vp)  =  -47i6(x-xI) 

where  k  is  the  wave  number,  p  is  density,  x  is  the  position  vector,  and  Xj  is  the 
spherical  point  source  position  vector.^  Semi-analytical  methods  such  as 
geometric  optics  (i.e.  ray  theory)  and  separation  of  variables  (i.e.  normal  modes  or 
fast  field),  as  well  as  pure  numerical  schemes  such  as  the  parabolic  equation 
methods,  can  be  used.  However,  as  diffraction  and  dispersion  can  become 
important  at  low  frequencies,  transfer  functions  obtained  using  full  wave  models, 
such  as  normal  mode  or  parabolic  equation  methods,  are  better.  Normal  mode 
theory  is  most  easily  used  in  ocean  environments  which  are  horizontally  stratified. 
Methods  exist  which  account  for  mode  coupling  effects  resulting  from  non¬ 
horizontal  stratification  and  rough  bathymetry.’^-*^  However,  transfer  functions 
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created  using  normal  mode  methods  may  require  significant  processing  time  and 
storage  for  deep  water  regions  where  many  modes  are  excited.  Additionally, 
separability  of  the  bottom  boundary  condition  can  become  questionable  if  rough 
bathymetry  exists. 

In  the  past,  the  parabolic  equation  method  has  been  limited  to  modelling 
narrow  angle  sound  energy  propagation.  However,  current  parabolic  acoustic 
models  have  been  able  to  model  sound  with  propagation  angles  close  to  90 
degrees,^*^^  Furthermore,  rough  bathymetry  and  variable  stratification  are  easily 
implemented.  Transfer  functions  constructed  using  parabolic  approximations  are 
especially  suited  for  low  frequency,  shallow  water,  rough  bathymetry  ocean 
conditions. 

3.  PARABOLIC  APPROXIMATION 

Parabolic  equation  methods  represent  the  acoustic  pressure  as  the  product 
of  two  functions: 

p(r,z)=  u(r,z)v(r)  (7) 

where  r  is  the  range  from  the  source  to  the  receiver,  so  thatr=rr-ro. 
Computation  time  for  a  numerical  solution  is  then  reduced  by  expressing  rapid 
variations  with  range  analytically,  and  leaving  the  slower  variations  in  range  to  be 
determined  numerically.  The  analytic  function  v(r)  changes  rapidly  with  range, 
while  the  numerical  function,  u(r,z),  slowly  varies  in  range  but  can  vary  rapidly 
with  depth. 

Substituting  (7)  into  the  Helmholtz  equation  and  neglecting  the  the 
backscattered  field  results  in  a  separable  equation  for  the  function  v(r)  which  has 
an  exponential  form.  FoUowing  Collins”,  v(r)  is: 
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(8) 


v(r)=;^i  kor 

Substituting  for  v(r)  and  factoring  results  in  a  parabolic  equation  for  u(r,z)  which 
includes  an  unwieldy  square  root  operator: 

^=iko(friX-  l)u 


where 


,2  i,2^  9^  9 


X  =  ^k^.lcd+  — -  f 


3z3z 


The  square  root  operator  has  been  simplified  using  either  a 
Taylor  series^^,  a  rational  function^*  or  a  Fade  series^'^^  approximation.  The  Finite 
Element  Parabolic  Equation  (FEPE)  of  Collins  uses  a  family  of  Fade  series  to 
create  higher  order  parabolic  equations  which  are  accurate  for  propagation 
angles  close  to  90  degrees.^'^^ 

Using  Fade  series,  the  FEPE  approximates  the  square  root  operator  as: 

iTTx-ut-^^ 

1  .1  v 


j=i 

a  n  = 

2n+l  2n+l 


(10) 


bj,n  - 


=  COS^ 

2n+l 


The  number  of  Fade  terms  used  determines  the  angle  of  sound  propagation  which 
can  be  accurately  modelled.  While  the  use  of  one  term  is  equivalent  to  the  wide 
angle  PE,  having  an  approximately  40  degree  propagation  angle,  the  use  of  four 
terms  provides  accuracy  of  close  to  90  degrees  depending  on  the  ocean 
environment.  The  amount  of  computer  time  required  increases  proportionally  with 
the  number  of  Fade  terms  used.^'^^ 

Substituting  the  parabolic  form  for  the  frequency  domain  transfer  function 
(Eq.  7),  the  time  domain  pressure  signal  can  be  evaluated  as: 
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p(rr,z,t)=  P(ro,zo,Ou(r,z,f)v(r,f)e-^"f‘df 


(11) 


Substituting  for  the  chosen  value  of  v(r) 


rr,z,t)=  J 


P(ro ,  zo ,  f)u(r,  z, 


Using  the  fact  that  ko  =  and  combining  exponential  terms,  (12)  becomes: 

CQ 


(12) 


p(rr ,  z,  t)  =  P(ro ,  Zo ,  f)u(r,  z,  f);^xp(-i27rf(t  -  -^j)  df  (13) 


With  a  change  of  time  variable: 

^0 

where  ^  is  the  reference  travel  time  of  a  plane  wave  solution,  the  resulting 


pressure  is 


rc 

rr,z,t')  = 


P(ro,  z,  f)u(r,  z,  f}^ 


(14) 


The  integral  expression  is  solved  numerically  using  the  inverse  fast  Fourier 
transform  (IFFT). 
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m.  MODEL  USE 


1.  MODEL  INPUT 

Modifications  in  FEPE_SYN,  a  prototype  Fourier  synthesis  code 
developed  by  Collins^* allowed  for  its  use  in  modelling  the  acoustic 
propagation  from  the  transmitter  to  Station  J  in  the  Monterey  Tomography 
Experiment.  The  Blackman  function  was  used  to  represent  the  transmitting  signal 
in  the  frequency  domain.  The  center  frequency  of  the  function  was  224  Hz,  the 
same  as  the  actual  source,  with  the  function  having  zero  values  below  212  Hz 
and  above  236  Hz  to  represent  the  source’s  ±  12  Hz  bandwidth  (Fig,  4).  The 
Blackman  function  provided  a  sharp  but  smooth  cutoff  above  and  below  the 
bandwidth,  reducing  the  sidelobes  in  the  time  domain.  The  inverse  Fourier 
transform  of  the  Blackman  function  closely  approximated  the  output  of  the 
correlator  receiver  used  in  the  experiment. 

The  bathymetry  used  was  obtained  from  NOAA^^  and  had  a  grid  spacing 
of  250  meters.  The  SVP  used  was  obtained  near  the  transmitter  during  the  1988 
Monterey  Bay  Experiment.  The  small  variation  between  the  sound  velocity 
profile  measured  at  the  source  and  the  receiver  allowed  the  deep  water  SVP  to  be 
used  along  the  entire  sound  path  (Fig.  5).  Bottom  sound  speed,  density  and 
attenuation  were  taken  to  be  the  average  values  for  continental  shelf  silty  sand  as 
obtained  by  Clay  and  Med  win 

The  sampling  interval  in  the  computational  example  is  one  millisecond  (ms). 
The  resulting  Nyquist  frequency  of  500  Hz,  is  more  than  twice  the  maximum 
frequency  used  of  236  Hz.  The  time  window  starts  440  ms  before  the  expected 
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arrival  of  a  plane  wave  solution,  travelling  at  a  speed  of  cq,  and  arriving  at  time 
t  =  r/cQ.  The  arrivals  are  tracked  for  1600  ms  after  the  reference  plane  wave 
solution  arrival. 

2.  COMPARISON  OF  MODEL  OUTPUT  WITH  DATA 

The  resulting  wide  angle,  one  Fade  term,  computer  modelled  output 
accounts  for  sound  energy  propagating  accurately  up  to  40  degrees,  and  shows  a 
series  of  arrivals  starting  140  ms  after  the  plane  wave  reference  (Fig.  6).  A  large 
double  peaked  pulse  is  shown  between  two  smaller  peaks.  Taking  the  higher 
order  terms  into  account,  using  two  Fade  terms,  the  constructive  and  destructive 
interference  of  the  higher  order  energy  is  seen  as  an  increase  in  the  amplitude  of 
the  predicted  signal,  and  a  narrowing  of  the  double  peaked  pulse  to  a  single  peak. 
The  major  pulse  amplitude  is  also  increased  (Fig.  7).  The  computer  output  from 
the  higher-order  example  shows  good  qualitative  agreement  with  the  measured 
station  J  data. 

The  predicted  and  the  measured  results  show  an  initial  arrival,  followed  by 
a  narrow,  larger  pulse  quickly  followed  by  a  series  of  smaller  amplitude  arrivals. 
The  spacing  between  the  pulses  and  the  relative  amplitudes  of  the  arrivals,  both 
plotted  as  amplitude  squared  in  Fig.  8,  are  in  agreement.  The  excellent  agreement 
of  the  higher-order  model  output  with  the  measured  data,  versus  that  of  the  lower 
order  output,  demonstrates  the  importance  of  including  the  higher-order  terms  in 
modelling  the  Monterey  Bay  propagation.  Higher-order  terms  can  be  important  in 
cases  where  steep  bathymetry  exists  such  as  in  this  example. 

Model  output  shows  a  sensitivity  to  the  horizontal  range  step  chosen.  The 
reproduction  of  the  measured  results,  by  the  time  domain  model,  required  the  use 
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Figure  4.  Frequency  domain  source  function  used  by  the  Finite  Equation 
Parabolic  Equation  (FEPE)  Fourier  synthesis  parabolic  equation  model  for 

the  Monterey  Bay  Experiment. 


Sound  Velocity  Profile  •  Source  lo  Station  J 


Figure  5.  Sound  velocity  profile  obtained  from  CTD  measurements  on 
December  14, 1988  near  the  source. 
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PE  Model  Output  tt  Station  I  -  15  tn  Range  Step 
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Figure  6.  Magnitude  squared  model  output  at  Station  J  using  a  15  m 
horizontal  range  step  and  one  Pad6  term. 
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Figure  7.  Magnitude  squared  model  output  at  Station  J  using  a  15  m  range 

step  and  two  Pad4  terms. 
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Figure  8.  Comparison  of  magnitude  squared  model  output  and  measured  data 

for  Station  J. 
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PR  Model  Ouipiii  at  Station  i  -  25  m  Range  Step.  2  Pad#  Tciins 
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Figure  9.  Magnitude  squared  model  output  at  Station  J  using  a  25  m  range 

step  and  two  Pad6  terms. 
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IV.  NORMAL  MODE  APPLICATIONS 


L  MODAL  DECOMPOSITION 

Modal  decomposition  provides  a  physically  meaningful  display  of  output 
from  the  parabolic  approximation.  In  addition,  normal  mode  travel  time 
information  can  be  used  for  tomographic  inversion.**^’^  In  particular,  time  domain 
modal  fluctuations  can  be  used  to  determine  the  surface  frequency  spectrum,  as 
verified  by  surface  buoys  deployed  in  the  1988  Monterey  Bay  Experiment. 

Sound  energy  arriving  at  a  receiver  can  be  described  as  the  sum  of  normal 
modes  (Fig.  10)  Zm(z): 

oo 

P(r,z,t)  =  X  am(r,t)  Zm(r,z,t)  (15) 

m=l 

where  p  is  pressure,  am  is  the  coefficient,  Zm  is  the  eigenfunction,  and  m  is  the 
mode  number.  Although  there  are  an  infinite  number  of  modes  possible,  only 
finite  number  of  modes  will  be  trapped  in  the  ocean  wave  guide,  depending  on 
the  frequency,  the  depth  of  the  channel,  the  SVP  and  the  boundary  conditions.^® 

The  number  of  trapped  modes  at  the  receiver  can  be  estimated  by  approximating 
the  coastal  shelf  area  of  Station  J  as  an  isospeed  sound  channel  bounded  by  a 
pressure  release  surface  and  a  sediment  bottom. 

Thirteen  trapped  modes  are  estimated  using  the  equation: 

m  <  ^2hcos  0c)  + 1  (16) 

where  m  is  the  number  of  modes  trapped,  f  is  frequency,  cq  is  the  sound  speed  in 
the  channel,  h  is  the  water  column  depth.^®  The  critical  angle,  0c,  is  defined  by 
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sin6c  _  1 
Co  Cl 


(17) 


where  ci  is  the  sediment  sound  speed. 

The  modes  are  computed  by  solving: 

(|j-k^-kS)z„,=0  (18) 

subject  to  the  boundary  conditions.  The  actual  SVP  and  a  finite  difference 
approximation  are  used. 

Since  Zm 's  are  orthogonal  functions,  i.e.. 


roo 


Zm(r,  z)Zn(r,  z)dz  =  6mn(z)  (19) 


JQ 

the  coefficient  am's  can  be  obtained  by  integrating  the  products  of  the  computed 
pressure  and  the  normal  modes: 


am=  I  P(n  ^  0-  Zm(r,  z)dz  (20) 

Jo 


The  pressure  time  series  can  be  seen  at  the  Station  J  receiver  for  each  mode  in  Fig. 

11. 


2.  MODEroENTIFICATION 

Acoustic  tomography  requires  that  measured  sound  arrivals  be  identified 
with  the  predicted  arrivals.  The  maximum  amplitude  measured  arrival  is  coincident 
with  the  maximum  predicted  sound  arrival.  The  sound  arrival  at  station  J  can  be 
seen  in  Fig.  1 1  to  be  clearly  dominated  by  mode  6  sound  energy.  The  second 
largest  amplitude  arrival,  which  just  precedes  the  largest  arrival,  is  also  dominated 
by  mode  six.  However,  modes  one,  five  and  eight  are  also  represented.  The  use  of 
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a  hydrophone  array  instead  of  a  single  receiver  should  make  it  possible  to 
measure  individual  modal  arrivals. 

3.  MODE  COUPLING 

Sound  energy  can  be  seen  to  move  from  higher  to  lower  modes  as  sound 
moves  up  the  slope  in  Monterey  Bay  to  Station  J  on  the  continental  shelf,  (Fig.  12 
and  Fig.  13).  Relating  normal  modes  to  rays  facilitates  the  explanation  of  this 
occurrence.  In  this  context,  ray  behavior  can  be  seen  as  the  constructive 
interference  of  neighboring  normal  modes.  Parameterizing  mode  quantities  by 
mode  number  and  equivalent  ray  angle  allows  an  "equivalent-ray"  to  be  found 
such  that  each  resonance  is  centered  at  an  equivalent-ray  angle.  These  rays  are 
"fuzzy"  in  that  they  have  finite  width,  and  then  a  sound  beam  may  be  used  to 
represent  a  ray  at  a  finite  frequency.^ The  ray  incident  angles  increase  as 
sound  travels  up  the  slope  due  to  bottom  interactions.  Ray  equivalent  angles  of 
modes  also  increase  as  the  depth  decreases. 

Mode  energy  is  concentrated  in  higher  modes  at  the  base  of  the  slope, 
with  mode  18  dominant.  The  equivalent-ray  angle  is  approximately  2.5  degrees, 
as  shown  in  Fig.  14.  A  ray  propagating  at  this  angle  will  interact  with  the  bottom 
only  once  on  the  slope.  The  bottom  interaction  will  increase  the  incident  angle  at 
the  continental  shelf  to  approximately  12  degrees.  The  mode  which  has  the 
resulting  ray-equivalent  angle  of  12  degrees  is  mode  six.  Therefore,  we  would 
expect  that  mode  6  energy  will  be  dominant  as  is  explained  using  a  ray  analogy. 
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Figure  10.  Eigenfunctions  calculated  for  Station  J.  The  sound  velocity  profile 

is  also  included. 


Figure  11.  The  magnitude  of  the  model  output  pressure  at  the  Station  J 
receiver  as  a  function  of  time  and  normal  mode  number. 


function  of  mode  number  at  ranges  45.75  and  49.25  km. 
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Figure  12c.  The  magnitude  squared  for  coefTicient  am  time  series  as  a  function 
of  mode  number  at  ranges  52.50  and  55.75  km  (Station  J). 
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Depth  (km) 


Figure  13.  Bathymetry  from  the  source  to  Station  J. 
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Angle  (deg) 


Range  (km) 


Figure  14.  Equivalent-ray  angles  as  a  function  of  mode  number  and  range. 

I 
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V.  CONCLUSIONS 


Fourier  synthesis  using  the  Collins- Westwood  parabolic  equation  model  is 
effective  in  modeling  short-duration  time  domain  signals.  Time  domain  parabolic 
equation  models  can  be  used  effectively  in  acoustically  complex  areas,  such  as 
Monterey  Canyon.  Modeling  propagation  in  areas  with  steep  slopes  require  a 
method  which  includes  higher  order  terms  and  which  conserves  energy. 

Model  output  from  the  Fourier  synthesis  parabolic  equation  model  may  be 
decomposed  into  normal  modes.  This  provides  a  tool  for  understanding  acoustic 
sound  transmission.  The  use  of  a  multiple  hydrophone  array  should  allow 
identification  of  individual  modal  sound  arrivals,  thereby  going  a  long  way 
toward  solving  the  forward  problem  for  acoustic  tomography  in  shallow  water. 
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